#include <btl/Maths/LapackSVD.hpp>

// declare the relevant LAPACK method
extern "C" {
   void dgesvd_(
         const char* JOBU, const char* JOBVT, int* M, int *N, double* A, int* lda,
         double* S, double *U, int* ldu, double* VT, int* ldvt,
         double* WORK, int* lwork, int* INFO);
}

namespace btl
{
namespace maths
{

void LapackSVD::doCompute(Eigen::MatrixXd& mat, int opts)
{
   _opts = opts;

   // need a full copy of the data
   double* const a = mat.data();
   int m = mat.rows();
   int n = mat.cols();

   _sv.resize(std::min(m,n), 1);
   _U.resize(m,m);
   _VT.resize(n,n);

   double* const sv = _sv.data();
   double* const u = _U.data();
   double* const vt = _VT.data();

   // ld* is the "leading" dimension (in other words, the stride)
   // of the relevant matrices (in column major this is the row count)
   int lda = m, ldu = m, ldvt = n;

   // input flags
   char JOBU = 'N', JOBVT = 'N';
   if ((opts & Eigen::ComputeThinU) == Eigen::ComputeThinU) JOBU = 'S';
   if ((opts & Eigen::ComputeFullU) == Eigen::ComputeFullU) JOBU = 'A';
   if ((opts & Eigen::ComputeThinV) == Eigen::ComputeThinV) JOBVT = 'S';
   if ((opts & Eigen::ComputeFullV) == Eigen::ComputeFullV) JOBVT = 'A';

   // output flags
   int LWORK;
   int INFO;

   // query LAPACK to find out how much working space it wants
   double size;
   LWORK = -1;
   ::dgesvd_(
      &JOBU, &JOBVT, &m, &n,
      a, &lda, sv, u, &ldu, vt, &ldvt,
      &size, &LWORK, &INFO);

   LWORK = (long int)(size);
   double* wk = new double[LWORK];

   ::dgesvd_(
      &JOBU, &JOBVT, &m, &n,
      a, &lda, sv, u, &ldu, vt, &ldvt,
      wk, &LWORK, &INFO);

   delete[] wk;

   // FIXME: check INFO here to see if there was a proble
}

} // namespace maths
} // namespace btl
